home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 188_01 / translib.for < prev    next >
Encoding:
Text File  |  1987-09-29  |  3.9 KB  |  161 lines
  1. /*
  2. HEADER:        ;
  3. TITLE:        C elementary transcendentals;
  4. VERSION:    1.0;
  5.  
  6. DESCRIPTION:       Source    code    for   all   standard    C                            
  7.                transcendentals
  8.  
  9.         Employs  ldexp()  and  frexp()  functions;  if                
  10.                suitable  versions  of these are  not  provided                
  11.                by  a given compiler,  the versions provided in                
  12.                source  code  will require  adaptation  to  the                
  13.                double float formats of the compiler.
  14.  
  15. KEYWORDS:    Transcendentals, library;
  16. SYSTEM:        CP/M-80, V3.1;
  17. FILENAME:    TRANS.C;
  18. WARNINGS:      frexp()   and  ldexp()  are   implementation                
  19.                dependent.   The  compiler  employed  does  not 
  20.                support    minus   (-)   unary   operators   in                
  21.                initializer  lists,  which are required by  the                
  22.                code.
  23. AUTHORS:    Tim Prince;
  24. COMPILERS:    MIX C, v2.0.1;
  25. */
  26.     program transt
  27.     real log
  28.     sum=0
  29.     do 1 i=1,1000
  30.       x=i*.01
  31. 1    sum=abs(tan(atan(x))/x-1)+sum
  32.     write(3,2)sum
  33. 2    format(' built-in tan(atan) error:',g15.7)
  34.     sum=0
  35.     do 3 i=1,1000
  36.       x=i*.01
  37. 3    sum=abs(exp(log(x))/x-1)+sum
  38.     write(3,4)sum
  39. 4    format(' built-in exp(log) error:',g15.7)
  40.     sum=0
  41.     do 5 i=1,1000
  42.       x=i*.01
  43. 5    sum=abs(xtan(xatan(x))/x-1)+sum
  44.     write(3,6)sum
  45. 6    format(' tan(atan) error:',g15.7)
  46.     sum=0
  47.     do 7 i=1,1000
  48.       x=i*.01
  49. 7    sum=abs(xexp(xlog(x))/x-1)+sum
  50.     write(3,8)sum
  51. 8    format(' exp(log) error:',g15.7)
  52.     end
  53.  
  54.       function xtan(x)
  55.     xtan=xsin(x)/xcos(x)
  56.     return
  57.       end
  58.  
  59.       function xcos(x)
  60.     xcos=xsin(1.57079633-x)
  61.     return
  62.       end
  63.  
  64.       function xsin(x)
  65.     real table(5)
  66.     data pi/3.141592653589793238/,
  67.      &      table/2.60198e-6,-.0001980746,.0083330258,
  68.      &        -.16666657,.99999999/
  69. c series to 28 bit accuracy
  70. c nearest multiple of pi
  71.     pim=anint(x/pi)
  72. c get remainder
  73.     xr=x-pi*pim
  74. c sin(x-pim*pi) = (-1)**pim * sin(x)
  75.     hpim=scal2(pim,-1)
  76. c if pim/2 not integer, pim is odd so change sign
  77.     if(aint(hpim).ne.hpim)xr=-xr
  78. c evaluate Horner polynomial, odd terms
  79.     xsin=poly(xr*xr,5,table)*xr
  80.     return
  81.       end
  82.  
  83.       function scal2(x,i)
  84. c    scal2=x*2**i
  85.     integer*1 ix(4),i
  86.     equivalence(xi,ix(1))
  87.     xi=x
  88.     ix(4)=i+ix(4)
  89.     scal2=xi
  90.     return
  91.       end
  92.  
  93.       function poly(x,n,table)
  94.     real table(n)
  95.     poly=table(1)
  96.     do 1 i=2,n
  97. 1    poly=poly*x+table(i)
  98.     return
  99.       end
  100.  
  101.       function xatan(x)
  102.     logical invert
  103.     real table(8)
  104.     data table/-.0046930,.0242506,-.0594827,.0991394,
  105.      &      -.1401932,.1996969,-.3333199,.9999999/
  106.     xr=abs(x)
  107.     invert=xr.gt.1
  108.     if(invert)xr=1/xr
  109.     xr=poly(xr*xr,8,table)*xr
  110. c tan(1/x)= tan(pi/2-x)
  111.     if(invert)xr=1.57079633-xr
  112.     xatan=sign(xr,x)
  113.     return
  114.       end
  115.  
  116.       function xlog(x)
  117.     xlog=alog2(x)*.693147181
  118.     return
  119.       end
  120.  
  121.       function alog2(x)
  122.     real table(3)
  123.     integer*1 ix(4),iexp
  124.     equivalence(xr,ix(1))
  125. c polynomial to 24 bit precision
  126.     data table/.5957779,.9615884,2.8853901/
  127.     xr=x
  128. c z'35' is leading 7 bits (minus hidden bit) of sqrt(.5)
  129. c shift so xr close to 1.
  130. c .true. value -1 on this system, bias z'80'
  131.     iexp=(ix(3).le.z'35')+ix(4)-z'80'
  132. c subtract 1 first so no accuracy lost
  133.     xr=scal2(xr,-iexp)-1
  134.     xr=xr/(xr+2)
  135. c polynomial in (xr-1)/(xr+1), .7 < xr < 1.4
  136.     alog2=poly(xr*xr,3,table)*xr+iexp
  137.     return
  138.       end
  139.  
  140.       function xexp(x)
  141.     xexp=exp2(x*1.442695040888963407)
  142.     return
  143.       end
  144.  
  145.       function exp2(x)
  146.     integer*1 expn
  147.     real table(7)
  148. c 26 bits precision -.5 < xr < .5 for 2**x
  149. c sinh (|x|<.7) use odd terms (x*1.442695)
  150.     data table/.000154,.00133904,.00961814,.05550358,
  151.      &      .2402265,.69314718,1./
  152.     if(x.ge.-128)go to 1
  153.       exp2=0
  154.       return
  155. 1    xr=amin1(x,126.9999)
  156.     expn=anint(xr)
  157.     xr=xr-expn
  158.     exp2=scal2(poly(xr,7,table),expn)
  159.     return
  160.       end
  161.